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IN PLANE STRESS ANALYSIS 


Summary 

In light of the discussions concerning the "blowing up" of the UK 
numerical results for the inplane stress calculations for silicon ribbon, 

I have decided to prepare this report that ties (I hope) the loose ends 
together. 

The results obtained since January 1985 show that the "blowing up" 
of computer solutions at least for smooth thermal Drofiles, is indeed 
a real physical phenomena. It is our view that they are precisely the 
"Luder's band" described in the Mobil Solar Quarterly report to JPL 
for January through March, 1984. 

ANALYSIS 

The inplane stresses that exist in a thin plate are governed by two 
general equations (equilibrium and compatability). Equilibrium is used in 
the form 
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which contains the normal stresses, a and a . The coordinate system used 

ax yy 

is shown in Fig. 1. Compatability is assumed in the form 
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which contains the (total) strain rates Physically compatability is 


h = N b e 
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when /J 2 > and the parameter J 2 is the second invariant of the deviatoric 

stress tensor. If D^N > /J 0 the value of h is zero. 

m Z 

The really novel and extremely important aspect of the Sumino-Haasen 
model is that the dislocation density changes where shown to vary according 
to the relationship. 
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We now postulate that this same law applies to the situation that occurs 
in ribbon growth. It will be shown below that having a changing dislocation 
density is a major effect in modeling real material behavior for the silicon 
used in solar cell applications. 

The major assumption that we are making here is that silicon is assumed 
to be isotropic in both its elastic and plastic constitutive relations. 

Eq. (6) happens to fit into the general framework of internal variable 

theories of viscoplasticity of Dillon and Kratochvil [3]. We sometimes change 
the parameter Kj in Eq. (6) to K^/10 in order to approximate a "shape factor" 

for the ribbon. We now assume thac the total strain rate c.. is the sum of 

the elastic and plastic parts. Thus our constitutive equation is 
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Steady State 

We henceforth assume that the growth process is In a steady state. Thus 
one can replace derivatives in time by those in the spatial coordinate x. 

By combining Eqs. (2) and (7), in the steady state we obtain 
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where v is the pull speed. In deriving Eq. (8) we neglected spatial deriva- 
tives of a and E. In the solution of Eqs. (8), we do allow a and E to vary 
with temperature and hence w'th x. In earlier numerical work where da/dx 
and dE/dx were retained this assumption was shown to be more than adequate. 

When the material is elastic, the solution of Eqs. (1) and (8) is readily 
accomplished on the digital computer for a specified thermal profile. We have 
also shown in previous reports that the more analytical procedures given in 
Boley and Weiner are not adequate for the thermal profiles of interest in the 
Mobil Solar and Westingnouse processes for growing silicon ribbon. In effect 

these profiles change too rapidly in space for the Boley-Weiner £4] method to 

2 

yield good results. On the other hand if V T is constant , the Boley-Weiner 
method gives excellent agreement with the stresses obtained using finite 
differences and the digital computer. The ribbon has no external loads applied 
to it so that stresses are due entirely to the thermal profile. 

Dislocation Density 

To this writer's knowledge this project is the first ever to use equa- 
tions such as Eq. (6) for any manufacturing process. That is we explicitly 
consider the changing internal structure and its effect on stresses and 
strains tnat develop during the growth of ribbon. 


Most materials researchers normally use materials with large values of 
the dislocation density. In this case changes in N do not drastically affect 
the material response. Some solar cell researchers however believe that 
electronic performance is reduced if the cells have a large number of 
active dislocations. Silicon people use words si.ch as "dislocation free" 
silicon for many applications. Certainly it appears that we should consider 
"low" dislocation densities and one is then forced to consider that they also 
change when silicon is subjected to stresses. In the growth processes, it is 
clear that a constitutive model with a changing dislocation density is an 

important improvement over a model that assumes N m remains constant. 

2 3 2 

Westinghouse claims average values of N are 10 to 10 /cm while Mobil 
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Solar material probably contains 10" to 10 /cm . These numbers apply to the 
respective final products in both cases. 

As will be discussed below, to study how dislocations grow one needs to 
assume the initial value of N at the melt interface , i.e. the value of N in 
the solid state at a temperature equal to the solidifying temperature. It 
turns out that this is a crucial assumption. We decided to try various values 
for at the melt interface and to evaluate the validity by comparing the 
predicted final values of N in the real product. Sometimes this is perfectly 
satisfactory approach, However we encountered "instabilities" in the numerical 
scheme and thus could not obtain valid solutions for all of the intended values 
of N q at the melt interface. This raised the question about the validity of our 
entire numerical code. 

Oxygen 

While we do not explicitly show how in Eqs. (4) and (6), our computer 
programs also allow for using different (but uniform) values for the oxygen 
content. Oxygen content is known to be different for float zone and CZ 
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which contains the normal stresses, a and a . The coordinate system used 
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which contains the (total) strain rates e^.. Physically compatability is 


silicon and to be a major actor in their having differing properties. It is 
probably different yet for the case of ribbon. 

The old EFG thermal profile used is 

T 4 = 437c' 1 - 36x cosnx + 1157e'- 066x 

- 317e"- 47 sin(^+ |) 

The modified EFG profile that is used is 

T 5 = 437e -1 ' 36 costtx + 1157e" ,066x 

- 158.5e"* 47x 


We continue to use 


T 3 = 1600e"* 0827x + 85e' 5x costtx 
+ 75 sin(TTX/4) - 273 

as the Westinghouse profile. 


RESULTS 

We assume thermal profiles which only vary with the distance x from the 
melt interface. We write Eqs. (1) and (8) in their finite difference approx- 
imations and treat the plastic incompatability integral in Eq. (8) as a pseudo 
thermal field. We then solve Eqs. (1) and (8) by an iterative procedure. The 



first Iteration Is obtained by assuming that the material is elastic, and 
hence the Integral In Eq. (8) is zero. The stresses so obtained are then 
used in Eqs. (5) and (6) to calculate the plastic strain rates and the dis- 
location density N m (x,y) as functions of spatial position. These new plastic 
strain rates are then used in the integral in Eq. (8), and the whole process 
repeated. In particular the a ■ stresses at x = 0 are determined in the 
standard way. They are very large and do not relax very much from their 
elastic values. 

Usually 20 iterations suffices for excellent results. Typically values 

O 

of /J 2 are of order 10 Pa, while changes in succeeding iterations are 
2 

usually only 10 Pa after 20 iterations. 

Typical results for the stresses and the dislocation density are shown 
in Figs. 2-4. The dislocation density changes by six (or more) orders of 
magnitude in Fig. (3). As shown in Fig. (2), the stresses eventually decrease 
in magnitude and N increases so much that of the ribbon is elastic and N m 
becomes zero. 

The data in Fig. (3) show that the value of the melt interface dislocation 

density that was used is (almost) ridicously small. On the other hand, the 

final value appears to be representative of ribbon products inferring at least 

that N q is perhaps not too bad. 

• PL 

The plastic strain rates e are shown in Fig. 4. Due to the low values 

of N q , the plastic strain rates are almost always small compared to the elastic 

strain rate. Moreover, one finds that the strain rates in the "cusp" region 

shown in Fig. 4 are very sensitive to the value of N Q assumed at the melt 

interface, as can be seen by comparing with the second curve given in Fig. 4. 

2 

It is relevant that the numerical scheme does not converge for N Q = 0.6/cm . 


A review of the numtri' J problem reveals that It is the plastic com- 
patability Integral In Eq. (8) which causes the solutions not to converge. 

This can be understood in view of the "cusp" region in Fig. 4 and noting that 
compatability requires spatial differentiation of the strain rates. When one 
considers the basic meaning of compatability, namely that the velocities are 
required to be single valued functions of space , one decides to look for 
experimental data for multivalued displacements or velocities. 

Such data exists in the Fig. 15 in the Mobil Solar report for January 
through March of 1984 £5] and has further been seen by us in EFG ribbon. In 
some cases one can actually see large "jumps" in the deformation field indica- 
ting multivalued displacements really do develop. According to personal 
communication with Dr. J. Kalejs of Mobil Solar, these regions of "Luder's 
bands" or some other forms of shear banding are developed at 1000°F at stresses 
of about 25 MPa. Our calculations, show the initial strain rate cusp develops 
at 1100°F or below and at stresses of 19.7 MPa or above. 

EFG Thermal Profile 

At this writing, we have obtained convergent solutions for the "old EFG" 
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thermal profile for the case where N q = 1,25 x 10 /cm . However we obtain 
nicely convergent solutions for a modified EFG profile for larger values of 
N q . The dislocation density at the melt interface for the modified EFG pro- 
file to result in convergent solutions is much below the value of N q used in 
the Westinghouse thermal profile. Numerical solutions for the EFG profile 
where N Q is close to the value which results in "unstable" solutions, shows 
two regions of "cusps" developing in plastic strain rate field in this thermal 
profile. The cusp near x = 0 is very large due to the very large values of 
a at x = 0 in these profiles. Considerably more care is involved in getting 

yy , 

stable solutions with an EFG thermal profile as compared to the Westinghouse 


one. 


Said differently, some computer solutions for EFG profile that blow-up 
are due to the spacing used. However this Is not normally the case with 
the Westinghouse profile. Very very small spacings are needed near x = 0 
to avoid blowing up for EFG. The root cause of this problem is the large 
stresses there. 

Inhomogeneojs plastic strain rat es 

Cusps in the plastic strain (rates) such as those shown in Fig. 4 occur 
also in the case of localized necking of ductile materials. Our localization 
is (probably) more accurately associated with shear binding than necking. 
However the analysis below, which parallels that used in necking, is believed 
to be useful in understanding the comp 1 exity of the phenomena shown in Fig. 4 
One solves Eq. (3) for the stress and obtains 
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The spatial derivative of the stress is 
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Since the process is a steady state one, we can write 3N/?x = N/y and 
• PL 

3t/3N, 3t/ 3T and 3 t/3e can be easily obtained from Eq. (9). Therefore we 
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have 



If one were to consider necking for example, In a constant temperature test 
In a material where the dislocation density Is uniform and constant force, 
then the R.H.S. of Eq. (11) Is zero and the plastic strain rate Is uniform 
In space as well . 

In the ribbon growing process however none of the terms on the R.H.S. 
of the analogue of Eq. (11) vanish and the strain rate in nonuniform is 
space. As shown in Fig. 4 there are regions where the strain rate is 
several orders of magnitude larger than it is elsewhere in the ribbon. There 
is no real simple interpretation of what causes the cusps to develop. Each 
term in the analogue of Eq. (11) makes a contribution that adds together 
to cause the "cusps". In the Westinghouse profile with a low melt face 
dislocation density, N is also small near x = 0 and the first and last terms 
on the R.H.S. of Eq. (11) tend to cancel one another. This keeps both 
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e and 3e /3x small. As the dislocation density grows N does as well and 
at x = 2.5 cm tha cusp has developed as shown In Fig. 4. 

The stresses in the (old) EFG profile are very large and nonuniform and 
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this in turn requires very small values of N in order to keep e from 
getting la». * and therefore Inhomogeneous. Small values of N exist if very 
small values of N Q a»*e assumed at the melt interface. 

The stresses at x = 0 

Hutchinson and Lambropoulos [6] used a finite element method to calculat 
the stresses in ribbon. They assumed a constitutive equation of the power la 
type and thereby assume (in effect) that the dislocation density did not 
change. They also assumed = 0 at the melt interface, x = 0. Subsequent 
work by Hutchinson justified this assumption. 


Our analysis shows that the o stresses at x ■ 0, are In fact very large 
and do not relax very much. An examination of the reason for this reveais 
t.nat It is because of the low values of N Q assumed. Both analyses are in fact 
"correct" and the question is then which is applicable to the growth process. 

W* believe that our analysis is mere applicable. We believe that the dis- 
location density is Indeed sma 1 ! and therefore one cannot assume it remains 
constant. 

Adaptive grids 

In recent years, it has become good computing practice to let the solution 
of a problem help distribute the grid points in order to concentrate more 
points in regions where rapid changes occur. We have incorporated that 
practice in many of our calculations. Caution is required though lest, the 
aspect ratio of the grid becomes unsatisfactory. Precautions such as using 
adaptive grids change the exact values where solutions blow up, but not enough 
to warrant discussion, at least for a Westinghouse profile. 
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